We need a bunch of these…
First jsonlite for processing the colours data pulled from the website
library(jsonlite)
Next, a bunch of data munging packages from the ‘tidyverse’
library(tidyr)
library(plyr) # rbind.fill useful for ragged data with missing entries
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(magrittr)
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
##
## extract
library(tibble)
library(stringr)
Next, the basic R spatial packages
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.3.1, PROJ 8.1.0
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:magrittr':
##
## extract
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
library(tmap)
library(tmaptools) # for geocode_OSM
tmap_mode("view") # for web maps
## tmap mode set to interactive viewing
Finally, a range of interpolation tools
library(spatstat) # for IDW
## Loading required package: spatstat.data
## Loading required package: spatstat.geom
## spatstat.geom 2.2-2
##
## Attaching package: 'spatstat.geom'
## The following objects are masked from 'package:raster':
##
## area, rotate, shift
## Loading required package: spatstat.core
## Loading required package: nlme
##
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
##
## getData
## The following object is masked from 'package:dplyr':
##
## collapse
## Loading required package: rpart
## spatstat.core 2.3-0
## Loading required package: spatstat.linnet
## spatstat.linnet 2.3-0
##
## spatstat 2.2-0 (nickname: 'That's not important right now')
## For an introduction to spatstat, type 'beginner'
library(maptools)
## Checking rgeos availability: TRUE
library(akima) # for triangulation
library(fields) # for splines
## Loading required package: spam
## Loading required package: dotCall64
## Loading required package: grid
## Spam version 2.7-0 (2021-06-25) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
##
## Attaching package: 'spam'
## The following objects are masked from 'package:base':
##
## backsolve, forwardsolve
## Loading required package: viridis
## Loading required package: viridisLite
## See https://github.com/NCAR/Fields for
## an extensive vignette, other supplements and source code
See the Dulux website for what this is all about
First I had a poke around on the website to figure out where the colour details were to be found
colour_groups <- c("blues", "browns", "greens", "greys", "oranges",
"purples", "reds", "whites-neutrals", "yellows")
base_url <- "https://www.dulux.co.nz/content/duluxnz/home/colour/all-colours.categorycolour.json/all-colours/"
The loop on the next slide
bind_rowscolours <- list()
for (i in 1:length(colour_groups)) {
colour_group <- colour_groups[i]
json_url <- str_c(base_url, colour_group)
json_file_name <- str_c(colour_group, ".json")
json <- fromJSON(json_url, flatten = TRUE)
# make a local copy (just for convenience)
write_json(json, json_file_name)
# get the colours information and add to list
the_colours <- rbind.fill(json$categoryColours$masterColour.colours)
colours[[i]] <- the_colours
Sys.sleep(0.5) # pause to not annoy the the server
}
df_colours <- bind_rows(colours)
write.csv(df_colours, "dulux-colours-raw.csv", row.names = FALSE)
head(df_colours)
## id red green blue lrv baseId name woodType coats
## 1 149253 205 210 206 67 vivid_white Pukaki Quarter None NA
## 2 149254 220 240 242 86 vivid_white Canoe Bay None NA
## 3 149255 226 240 245 87 vivid_white Mt Dobson None NA
## 4 149256 220 230 235 80 vivid_white Raetihi None NA
## 5 149257 217 219 223 74 vivid_white Taiaroa Head None NA
## 6 149258 180 200 219 60 vivid_white Gulf Harbour None NA
There are paint names with modifiers as suffixes for different shades of particular colours, and we need to handle this
The modifiers are
paint_modifiers <- c("Half", "Quarter", "Double")
Here’s one way to clean this up (there are others…)
df_colours_tidied <- df_colours %>%
## remove some columns we won't be needing
select(-id, -baseId, -woodType, -coats) %>%
## separate the name components, filling from the left with NAs if <5
separate(name, into = c("p1", "p2", "p3", "p4", "p5"), sep = " ",
remove = FALSE, fill = "left") %>%
## replace any NAs with an empty string
mutate(p1 = str_replace_na(p1, ""),
p2 = str_replace_na(p2, ""),
p3 = str_replace_na(p3, ""),
p4 = str_replace_na(p4, "")) %>%
## if p5 is a paint modifiers, then recompose name
## from p1:p4 else from p1:p5
## similarly keep modifier where it exists
mutate(placename = if_else(p5 %in% paint_modifiers,
str_trim(str_c(p1, p2, p3, p4, sep = " ")),
str_trim(str_c(p1, p2, p3, p4, p5, sep = " "))),
modifier = if_else(p5 %in% paint_modifiers,
p5, "")) %>%
## remove some places that are tricky to deal with later
filter(!placename %in% c("Chatham Islands",
"Passage Rock",
"Auckland Islands",
"Cossack Rock")) %>%
## throw away variables we no longer and reorder
select(name, placename, modifier, red, green, blue)
# save it so we have it for later
write.csv(df_colours_tidied, "dulux-colours.csv", row.names = FALSE)
Add x and y columns to our data for the coordinates
df_colours_tidied_xy <- df_colours_tidied %>%
mutate(x = 0, y = 0)
tmaptools::geocode_OSMCode on the next slide
x y coordinates as we have space for (due to the modifiers) from the geocoding resultsBest not to re-run this (it takes a good 10 minutes and it’s not good to repeatedly geocode and hit the OSM server)
for (placename in unique(df_colours_tidied_xy$placename)) {
address <- str_c(placename, "New Zealand", sep = ", ")
geocode <- geocode_OSM(address, as.data.frame = TRUE, return.first.only = FALSE)
num_geocodes <- nrow(geocode)
matching_rows <- which(df_colours_tidied_xy$placename == placename)
for (i in 1:length(matching_rows)) {
if (!is.null(geocode)) {
if (num_geocodes >= i) {
df_colours_tidied_xy[matching_rows[i], ]$x <- geocode$lon[i]
df_colours_tidied_xy[matching_rows[i], ]$y <- geocode$lat[i]
}
}
}
Sys.sleep(0.5) # so as not to over-tax the geocoder
}
# Remove any we missed
df_colours_tidied_xy <- df_colours_tidied_xy %>%
filter(x != 0 & y != 0)
write.csv(df_colours_tidied_xy, "dulux-colours-xy.csv", row.names = FALSE)
Make the dataframe into a sf point dataset
dulux_colours_sf <- st_as_sf(df_colours_tidied_xy,
coords = c("x", "y"), # columns with the coordinates
crs = 4326) %>% # EPSG:4326 for lng-lat
st_transform(2193) %>% # convert to NZTM
## and make an RGB column
mutate(rgb = rgb(red / 255, green / 255, blue/ 255))
# jitter any duplicate locations
duplicate_pts <- which(duplicated(dulux_colours_sf$geometry) |
duplicated(dulux_colours_sf$geometry, fromLast = TRUE))
jittered_pts <- dulux_colours_sf %>%
slice(duplicate_pts) %>%
st_jitter(50)
dulux_colours_sf[duplicate_pts, ]$geometry <- jittered_pts$geometry
st_write(dulux_colours_sf, "dulux-colours-pts.gpkg", delete_dsn = TRUE)
Remember to update the dataframe of points with the jittered points
jittered_pts <- dulux_colours_sf %>%
st_coordinates() %>%
as_tibble()
df_colours_tidied_xy <- df_colours_tidied_xy %>%
mutate(x = jittered_pts$X, y = jittered_pts$Y)
write.csv(df_colours_tidied_xy, "dulux-colours-xy-jit.csv", row.names = TRUE)
nz <- st_read("nz.gpkg")
## Reading layer `nz' from data source
## `/home/osullid3/Documents/code/dulux-colours-map/nz.gpkg' using driver `GPKG'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1089972 ymin: 4748123 xmax: 2091863 ymax: 6223160
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
## Reading layer `dulux-colours-pts' from data source
## `/home/osullid3/Documents/code/dulux-colours-map/dulux-colours-pts.gpkg'
## using driver `GPKG'
## Simple feature collection with 904 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 1094875 ymin: 4759975 xmax: 2073657 ymax: 6191716
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
tm_shape(nz) +
tm_borders() +
tm_shape(dulux_colours_sf) +
tm_dots(col = "rgb") +
tm_basemap("Esri.WorldTopoMap")
Points aren’t really much fun, instead:
We can make Voronois and clip to NZ
dulux_colours_vor <- dulux_colours_sf %>%
st_union() %>%
st_voronoi() %>%
st_cast() %>%
st_as_sf() %>%
st_join(dulux_colours_sf, left = FALSE) %>%
st_intersection(st_read("nz.gpkg"))
st_write(dulux_colours_vor, "dulux-colours-vor.gpkg", delete_dsn = TRUE)
## Reading layer `dulux-colours-vor' from data source
## `/home/osullid3/Documents/code/dulux-colours-map/dulux-colours-vor.gpkg'
## using driver `GPKG'
## Simple feature collection with 902 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1089972 ymin: 4748123 xmax: 2091863 ymax: 6223160
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000
tm_shape(dulux_colours_vor) +
tm_polygons(col = "rgb", id = "placename",
alpha = 0.75, border.col = "grey", lwd = 0.2) +
tm_basemap("Esri.WorldTopoMap")
Using the akima::interp package
components = c("red", "green", "blue")
layers = list()
for (component in components) {
# the dimensions, nx, ny give ~2500m resolution
layers[[component]] <- raster(
interp(df_colours_tidied_xy$x, df_colours_tidied_xy$y,
df_colours_tidied_xy[[component]],
nx = 402, ny = 591, linear = TRUE))
}
rgb.t <- brick(layers)
crs(rgb.t) <- st_crs(nz)$wkt
rgb.t <- mask(rgb.t, nz)
writeRaster(rgb.t, "dulux-colours-tri.tif", overwrite = TRUE)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(rgb.t) +
tm_rgb() +
tm_basemap("Esri.WorldTopoMap")
Here we use the spatstat::idw function
We have to make spatstat::ppp (point pattern) objects
## We need a window for the ppps
W <- nz %>%
st_geometry() %>%
st_buffer(1000) %>%
st_union() %>%
as("Spatial") %>%
as.owin()
layers <- list()
for (component in components) {
pp <- ppp(x = df_colours_tidied_xy$x,
y = df_colours_tidied_xy$y,
window = W,
marks = df_colours_tidied_xy[[component]])
## eps is the approximate resolution
layers[[component]] <- raster(idw(pp, eps = 2500, power = 4))
}
rgb.idw <- brick(layers)
crs(rgb.idw) <- st_crs(nz)$wkt
writeRaster(rgb.idw, "dulux-colours-idw.tif", overwrite = TRUE)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(rgb.idw) +
tm_rgb() +
tm_basemap("Esri.WorldTopoMap")
For this we use the fields::Tps (thin plate spline) function
We need an empty raster to interpolate into
r <- nz %>%
st_make_grid(cellsize = 2500, what = "centers") %>%
st_coordinates() %>%
as_tibble() %>%
rename(x = X, y = Y) %>%
mutate(z = 0) %>%
rasterFromXYZ()
Then we can interpolate as before
layers <- list()
for (component in components) {
spline <- Tps(df_colours_tidied_xy[, c("x", "y")],
df_colours_tidied_xy[[component]],
scale.type = "unscaled", m = 3)
layers[[component]] <- interpolate(r, spline)
}
rgb.s <- brick(layers)
crs(rgb.s) <- st_crs(nz)$wkt
rgb.s <- mask(rgb.s, nz)
writeRaster(rgb.s, "dulux-colours-spline.tif", overwrite = TRUE)
## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tm_shape(rgb.s) +
tm_rgb() +
tm_basemap("Esri.WorldTopoMap")
## Warning: Raster values found that are outside the range [0, 255]
sf and tmap (for basic geospatial)